home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / libgutil / fft.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  5KB  |  301 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. /*
  18.  *    fft -
  19.  *        simple 2-D fft support
  20.  *
  21.  *            Paul Haeberli - 1990
  22.  */
  23. #include "math.h"
  24. #include "image.h"
  25. #include "fft.h"
  26.  
  27. imgtocomp(image,comp,size)
  28. IMAGE *image;
  29. complex *comp;
  30. int size;
  31. {
  32.     int x, y;
  33.     short *row;
  34.  
  35.     row = (short *)malloc(size*sizeof(short));
  36.     for(y=0; y<size; y++) {
  37.     getrow(image,row,y,0);
  38.     for(x=0; x<size; x++) {
  39.         comp->r = row[x]/255.0;
  40.         comp->i = 0.0;
  41.         comp++;
  42.     }
  43.     }
  44.     free(row);
  45. }
  46.  
  47. zapzero(comp,size)
  48. complex *comp;
  49. int size;
  50. {
  51.     int x, y;
  52.  
  53.     for(y=0; y<size; y++) {
  54.     for(x=0; x<size; x++) {
  55.         if(x == 0 || y == 0) {
  56.         comp->r = 0.0;
  57.         comp->i = 0.0;
  58.         }
  59.         comp++;
  60.     }
  61.     }
  62. }
  63.  
  64. #define WRAPX(x)    (((x)+size/2)%size)
  65. #define WRAPY(y)    (((y)+size/2)%size)
  66.  
  67. ffttoimg(buf,image,size)
  68. complex *buf;
  69. IMAGE *image;
  70. int size;
  71. {
  72.     complex *comp;
  73.     int x, y;
  74.     float min, max, sc, val;
  75.     int n, i;
  76.     float mag, *mags;
  77.     short *row;
  78.  
  79.     n = size*size;
  80.     mags = (float *)malloc(n*sizeof(float));
  81.     row = (short *)malloc(size*sizeof(short));
  82.     max = -100000.0;
  83.     min = 100000.0;
  84.     comp = buf; 
  85.     comp->r = 0;
  86.     comp->i = 0;
  87.     i = 0;
  88.     while(n--) {
  89.     mag = sqrt(comp->r*comp->r+comp->i*comp->i);
  90.     mags[i++] = mag;
  91.     if(max<mag)
  92.         max = mag;
  93.     if(min>mag)
  94.         min = mag;
  95.     comp++;
  96.     }
  97.     if(min<0.0)
  98.     min = -min;
  99.     if(max<0.0)
  100.     max = -max;
  101.     if(min>max)
  102.        max = min;
  103.     printf("min max %f %f\n",min,max);
  104.     sc = 1.0/max;
  105.     comp = buf; 
  106.     i = 0;
  107.     for(y=0; y<size; y++) {
  108.     for(x=0; x<size; x++) {
  109.         mag = mags[i++];;
  110.         row[WRAPX(x)] = 255.0*(mag-min)/(max-min);
  111.         comp++;
  112.     }
  113.     putrow(image,row,WRAPY(y),0);
  114.     }
  115.     free(row);
  116.     free(mags);
  117. }
  118.  
  119. complex *compalloc(n)
  120. int n;
  121. {
  122.     return (complex *)malloc(n*sizeof(complex));
  123. }
  124.  
  125. fftrange(comp,n,min,max)
  126. complex *comp;
  127. int n;
  128. double *min, *max;
  129. {
  130.     *max = -100000.0;
  131.     *min = 100000.0;
  132.     while(n--) {
  133.     if(comp->r > *max)
  134.         *max = comp->r;
  135.     if(comp->r < *min)
  136.         *min = comp->r;
  137.     comp++;
  138.     }
  139. }
  140.  
  141. comptoimg(comp,image,size)
  142. complex *comp;
  143. IMAGE *image;
  144. int size;
  145. {
  146.     int x, y;
  147.     double min, max;
  148.     short *row;
  149.  
  150.     row = (short *)malloc(size*sizeof(short));
  151.     fftrange(comp,size*size,&min,&max);
  152.     printf("final img range: %f %f\n",min,max);
  153.     for(y=0; y<size; y++) {
  154.     for(x=0; x<size; x++) {
  155.         row[x] = 255.0*comp->r;
  156.         if(row[x]>255)
  157.             row[x]=255;
  158.         if(row[x]<0)
  159.             row[x]=0;
  160.         comp++;
  161.     }
  162.     putrow(image,row,y,0);
  163.     }
  164.     free(row);
  165. }
  166.  
  167. printcomp(comp,size)
  168. complex *comp;
  169. int size;
  170. {
  171.     int x, y;
  172.  
  173.     for(y=0; y<size; y++) {
  174.     for(x=0; x<size; x++) {
  175.         if(comp->r != 0.0 || comp->i != 0.0) 
  176.         printf("at %d %d: %f %f\n",x,y,comp->r,comp->i);
  177.         comp++;
  178.     }
  179.     }
  180. }
  181.  
  182. fft2(n,isign,din,dout) 
  183. int n, isign;
  184. complex *din;
  185. complex *dout;
  186. {
  187.     int y;
  188.  
  189.     for(y=0; y<n; y++) 
  190.     fft(n,isign,din+(y*n));
  191.     ffttranspose(n,din,dout);
  192.     for(y=0; y<n; y++) 
  193.     fft(n,isign,dout+(y*n));
  194. }
  195.  
  196. ffttranspose(n,din,dout)
  197. int n;
  198. complex *din;
  199. complex *dout;
  200. {
  201.     int x, y;
  202.     complex *sptr, *dptr;
  203.  
  204.     sptr = din;
  205.     for(y=0; y<n; y++) {
  206.     dptr = dout+y;
  207.     for(x=0; x<n; x++) {
  208.         *dptr = *sptr++;
  209.         dptr += n;
  210.     }
  211.     }
  212. }
  213.  
  214. fft(n,isign,data)
  215. int n, isign;
  216. complex data[];
  217. {
  218.  
  219.     int i;
  220.     int j, k, ln, nv2;
  221.     complex temp, u, w;
  222.     int le, le1,l;
  223.  
  224.     powercheck(n);
  225.     if(isign>0)
  226.     isign = 1;
  227.     else
  228.     isign = -1;
  229.     for(ln=0, i=n; i>1; i>>=1)
  230.     ln++;
  231.     nv2 = n/2;
  232.     j=0;
  233.     n--;
  234.     for(i=0; i<n; i++) {
  235.     if(i<j) {
  236.         temp = data[j];
  237.         data[j] = data[i];
  238.         data[i] = temp;
  239.     }
  240.     k = nv2;
  241.     while(j>=k) {
  242.         j -= k;
  243.         k >>= 1;
  244.     }
  245.     j += k;
  246.     }
  247.     n++;
  248.     le = 2;
  249.     le1 = 1;
  250.     for(l=0; l<ln; l++) {
  251.     double tempr, tempi;
  252.     double ur, ui;
  253.     double wr, wi;
  254.     complex *ipptr, *iptr;
  255.  
  256.     ur = 1.0;
  257.     ui = 0.0;
  258.     wr = cos(isign*3.1415926535/le1);
  259.     wi = -sin(isign*3.1415926535/le1);
  260.     for(j=0; j<le1; j++) {
  261.         for(i=j; i<n; i+=le) {
  262.         iptr = data+i;
  263.         ipptr = iptr+le1;
  264.         tempr = ur*ipptr->r - ui*ipptr->i;
  265.         tempi = ur*ipptr->i + ui*ipptr->r;
  266.         ipptr->r = iptr->r - tempr;
  267.         ipptr->i = iptr->i - tempi;
  268.         iptr->r += tempr;
  269.         iptr->i += tempi;
  270.         }
  271.         tempr = ur;
  272.         ur = wr*ur - wi*ui;
  273.         ui = wr*ui + wi*tempr;
  274.     }
  275.     le <<= 1;
  276.     le1 <<= 1;
  277.     }
  278.     if(isign == 1) {
  279.     for(i=0; i<n; i++) {
  280.         data[i].r /= n;
  281.         data[i].i /= n;
  282.     }
  283.     } 
  284. }
  285.  
  286. powercheck(n)
  287. int n;
  288. {
  289.     int i;
  290.  
  291.     for(i=0; i<32; i++) {
  292.     if(n&1)
  293.         break;
  294.     n >>= 1;
  295.     }
  296.     if(n != 1) {
  297.     fprintf(stderr,"fft: n points must be a power of 2!!\n");
  298.     exit(1);
  299.     }
  300. }
  301.